import numpy as np
from scipy.spatial.distance import cdist
from sklearn.neighbors import KernelDensity
from collections import deque
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
import platform
import warnings

# ================== 跨平台字体配置 ==================
def config_visualization():
    """配置可视化参数（中文+数学符号）"""
    system = platform.system()
    font_config = {
        'Windows': {'zh': 'simhei.ttf', 'math': 'STIXGeneral'},
        'Darwin': {'zh': 'STHeiti Medium.ttc', 'math': 'Apple Symbols'},
        'Linux': {'zh': 'wqy-zenhei.ttc', 'math': 'DejaVu Sans'}
    }
    
    try:
        # 获取系统字体路径
        if system == 'Windows':
            zh_path = f"C:/Windows/Fonts/{font_config[system]['zh']}"
        elif system == 'Darwin':
            zh_path = f"/System/Library/Fonts/{font_config[system]['zh']}"
        else:
            zh_path = f"/usr/share/fonts/truetype/wqy/{font_config[system]['zh']}"

        # 注册中文字体
        zh_font = fm.FontProperties(fname=zh_path)
        plt.rcParams['font.family'] = zh_font.get_name()
        
        # 配置数学符号
        plt.rcParams['mathtext.fontset'] = font_config[system]['math']
        plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示
        
        # 设置备用字体栈
        plt.rcParams['font.sans-serif'] = [zh_font.get_name(), 'SimHei', 
                                         font_config[system]['math']]
        return zh_font
    
    except Exception as e:
        warnings.warn(f"字体配置异常: {str(e)}，使用默认配置")
        plt.rcParams['axes.unicode_minus'] = False
        return fm.FontProperties()

# 初始化字体配置
zh_font = config_visualization()

# ================== 混沌生成器（数值安全版）==================
def tent_chebyshev_map(n, dim, steps=100):
    """生成混沌序列（带数值安全限制）"""
    seq = np.zeros((steps, dim))
    x = np.random.uniform(0.2, 0.8, dim)  # 初始值避开敏感区域
    
    for i in range(steps):
        # Tent映射部分
        mask = x < 0.5
        x[mask] = (2 * x[mask]).clip(0.1, 0.9)
        x[~mask] = (2 * (1 - x[~mask])).clip(0.1, 0.9)
        
        # 添加随机扰动（限制幅度）
        delta = np.random.uniform(-0.005, 0.005, dim)
        x = (x + delta).clip(0.05, 0.95)
        
        # Chebyshev变换（带数值保护）
        x_scaled = 1.999 * x - 0.999  # 映射到[-0.999, 0.999]
        x_cheb = np.cos(3 * np.arccos(x_scaled.clip(-1,1)))
        x = (x_cheb + 1) / 2.0  # 重新归一化到[0,1]
        
        seq[i] = x
    
    return seq[-n:]

# ================== 核心算法类 ==================
class ECPSO:
    def __init__(self, obj_func, dim, bounds, pop_size=50, max_iter=500):
        # 初始化参数
        self.obj_func = obj_func
        self.dim = dim
        self.bounds = np.array(bounds)
        self.range = self.bounds[:,1] - self.bounds[:,0]
        self.pop_size = pop_size
        self.max_iter = max_iter
        
        # PSO参数
        self.omega = 0.729
        self.c1 = self.c2 = 1.494
        self.alpha0 = 0.1
        self.beta0 = 0.05
        self.window_size = 5
        
        # 初始化种群
        self.pop_pos = self._init_population()
        self.pop_vel = np.random.uniform(-1, 1, (pop_size, dim)) * 0.1 * self.range
        self.pbest_pos = self.pop_pos.copy()
        self.pbest_val = np.array([obj_func(x) for x in self.pop_pos])
        self.gbest_pos = self.pbest_pos[np.argmin(self.pbest_val)]
        self.gbest_val = np.min(self.pbest_val)
        
        # 历史记录
        self.history = deque(maxlen=10)
        self.entropy_window = deque(maxlen=self.window_size)
        
        # 熵计算器
        self.kde = KernelDensity(bandwidth=0.15)

    def _init_population(self):
        """安全初始化种群位置"""
        chaos = tent_chebyshev_map(self.pop_size, self.dim)
        return chaos * self.range + self.bounds[:,0]

    def optimize(self):
        """执行优化主循环"""
        convergence = []
        for iter in range(self.max_iter):
            # 速度更新（含历史动量）
            self._update_velocity(iter)
            
            # 位置更新与截断
            self.pop_pos = np.clip(self.pop_pos + self.pop_vel, 
                                 self.bounds[:,0], self.bounds[:,1])
            
            # 更新最优解
            self._update_bests()
            
            # 每5代执行混沌扰动
            if iter % 5 == 0:
                self._chaos_perturbation(iter)
                self._calculate_entropy()
            
            convergence.append(self.gbest_val)
        
        return self.gbest_pos, self.gbest_val, convergence

    def _update_velocity(self, iter):
        """更新粒子速度"""
        history_term = sum(0.8**k * (h - self.pop_pos) 
                        for k, h in enumerate(self.history))
        
        # 计算梯度场
        gradient = self._calculate_gradient()
        
        # 梯度敏感的混沌步长
        chaos_vel = tent_chebyshev_map(self.pop_size, self.dim)
        chaos_term = chaos_vel / (np.linalg.norm(gradient) + 1e-8)
        
        self.pop_vel = (self.omega * self.pop_vel +
                       self.c1 * np.random.rand() * (self.pbest_pos - self.pop_pos) +
                       self.c2 * np.random.rand() * (self.gbest_pos - self.pop_pos) +
                       0.3 * history_term +
                       0.1 * chaos_term)

    def _calculate_gradient(self):
        """计算目标函数的梯度"""
        gradient = np.zeros(self.dim)
        for d in range(self.dim):
            # 使用中心差分法计算梯度
            delta = np.zeros(self.dim)
            delta[d] = 1e-5 * self.range[d]
            pos_plus = self.gbest_pos + delta
            pos_minus = self.gbest_pos - delta
            gradient[d] = (self.obj_func(pos_plus) - self.obj_func(pos_minus)) / (2 * delta[d])
        return gradient

    def _update_bests(self):
        """更新个体和全局最优"""
        current_val = np.array([self.obj_func(x) for x in self.pop_pos])
        improved = current_val < self.pbest_val
        
        self.pbest_pos[improved] = self.pop_pos[improved]
        self.pbest_val[improved] = current_val[improved]
        
        if (min_val := np.min(current_val)) < self.gbest_val:
            self.gbest_val = min_val
            self.gbest_pos = self.pop_pos[np.argmin(current_val)]
            self.history.append(self.gbest_pos.copy())

    def _chaos_perturbation(self, iter):
        """应用混沌扰动"""
        alpha = self.alpha0 * np.exp(-iter / (0.2*self.max_iter))
        chaos = tent_chebyshev_map(3, self.dim) * 2 - 1  # [-1,1]范围
        new_pos = self.gbest_pos + alpha * chaos * self.range
        
        # 替换最差粒子
        new_val = np.array([self.obj_func(x) for x in new_pos])
        worst = np.argmax(self.pbest_val)
        if new_val.min() < self.pbest_val[worst]:
            self.pop_pos[worst] = new_pos[np.argmin(new_val)]
            self.pbest_val[worst] = new_val.min()

    def _calculate_entropy(self):
        """计算种群分布熵"""
        X_norm = (self.pop_pos - self.bounds[:,0]) / self.range
        self.kde.fit(X_norm)
        log_prob = self.kde.score_samples(X_norm)
        prob = np.exp(log_prob)
        entropy = -np.mean(prob * log_prob)
        self.entropy_window.append(entropy)
        
        # 动态窗口熵
        sliding_H = self._calculate_sliding_entropy()
        # 根据滑动窗口熵调整参数
        self.omega = 0.5 + 0.5 * (1 - sliding_H / self.max_entropy)
    
    def _calculate_sliding_entropy(self):
        """计算动态窗口熵"""
        window_entropies = list(self.entropy_window)
        if not window_entropies:
            return 0.0
        weights = np.exp(-np.arange(len(window_entropies))/5)
        weights /= weights.sum()
        sliding_H = np.dot(window_entropies, weights)
        self.max_entropy = max(self.max_entropy, sliding_H) if hasattr(self, 'max_entropy') else sliding_H
        return sliding_H

# ================== 动态环境模拟器 ==================
class DynamicPeaks:
    def __init__(self, dim, num_peaks=3, speed=0.1, bounds=(-5,5)):
        self.dim = dim
        self.bounds = np.array(bounds)
        self.peaks = [self._generate_peak() for _ in range(num_peaks)]
        self.velocities = [np.random.uniform(-speed, speed, dim) 
                          for _ in range(num_peaks)]
        self.current = 0
    
    def _generate_peak(self):
        """生成安全峰位置"""
        return np.random.uniform(self.bounds[0]+0.5, self.bounds[1]-0.5, self.dim)
    
    def shift(self):
        """移动峰位置"""
        for i in range(len(self.peaks)):
            self.peaks[i] = np.clip(self.peaks[i] + self.velocities[i],
                                   self.bounds[0], self.bounds[1])
        self.current = (self.current + 1) % len(self.peaks)
    
    def __call__(self, x):
        return np.sum((x - self.peaks[self.current]) ** 2)

# ================== 算法基类 ==================
class Optimizer:
    """优化算法基类"""
    def __init__(self, obj_func, dim, bounds, pop_size=50, max_iter=500):
        self.obj_func = obj_func
        self.dim = dim
        self.bounds = np.array(bounds)
        self.range = self.bounds[:,1] - self.bounds[:,0]
        self.pop_size = pop_size
        self.max_iter = max_iter
        self.pop_pos = None
        self.gbest_val = float('inf')
        self.gbest_pos = None

    def _init_population(self):
        """初始化种群"""
        self.pop_pos = np.random.uniform(self.bounds[:,0], self.bounds[:,1], 
                                       (self.pop_size, self.dim))
    
    def optimize(self):
        """优化流程模板"""
        raise NotImplementedError

# ================== 标准PSO ==================
class PSO(Optimizer):
    """标准粒子群算法"""
    def __init__(self, obj_func, dim, bounds, pop_size=50, max_iter=500):
        super().__init__(obj_func, dim, bounds, pop_size, max_iter)
        self.omega = 0.729
        self.c1 = self.c2 = 1.494
        self.pbest_pos = None
        self.pbest_val = None
        self.vel = None
        
    def _init_velocity(self):
        """初始化速度"""
        self.vel = np.random.uniform(-1, 1, 
                                   (self.pop_size, self.dim)) * 0.1 * self.range

    def optimize(self):
        self._init_population()
        self._init_velocity()
        
        # 初始化最优解
        self.pbest_pos = self.pop_pos.copy()
        self.pbest_val = np.array([self.obj_func(x) for x in self.pop_pos])
        self.gbest_val = np.min(self.pbest_val)
        self.gbest_pos = self.pbest_pos[np.argmin(self.pbest_val)]
        
        convergence = []
        for _ in range(self.max_iter):
            # 更新速度
            r1, r2 = np.random.rand(2)
            self.vel = (self.omega * self.vel +
                       self.c1 * r1 * (self.pbest_pos - self.pop_pos) +
                       self.c2 * r2 * (self.gbest_pos - self.pop_pos))
            
            # 更新位置
            self.pop_pos = np.clip(self.pop_pos + self.vel, 
                                 self.bounds[:,0], self.bounds[:,1])
            
            # 评估
            current_val = np.array([self.obj_func(x) for x in self.pop_pos])
            
            # 更新最优
            improved = current_val < self.pbest_val
            self.pbest_pos[improved] = self.pop_pos[improved]
            self.pbest_val[improved] = current_val[improved]
            
            if np.min(current_val) < self.gbest_val:
                self.gbest_val = np.min(current_val)
                self.gbest_pos = self.pop_pos[np.argmin(current_val)]
            
            convergence.append(self.gbest_val)
        
        return self.gbest_pos, self.gbest_val, convergence

# ================== 混沌PSO ================== 
class CPSO(PSO):
    """混沌粒子群算法"""
    def __init__(self, obj_func, dim, bounds, pop_size=50, max_iter=500):
        super().__init__(obj_func, dim, bounds, pop_size, max_iter)
        self.chaos_steps = 5
        
    def _chaos_perturb(self):
        """混沌扰动"""
        chaos = tent_chebyshev_map(self.pop_size, self.dim, self.chaos_steps)
        chaos = chaos * self.range + self.bounds[:,0]
        chaos_val = np.array([self.obj_func(x) for x in chaos])
        
        # 替换最差粒子
        worst = np.argmax(self.pbest_val)
        if chaos_val.min() < self.pbest_val[worst]:
            self.pop_pos[worst] = chaos[np.argmin(chaos_val)]
            self.pbest_val[worst] = chaos_val.min()

    def optimize(self):
        convergence = super().optimize()[2]
        for iter in range(self.max_iter):
            if iter % 10 == 0:
                self._chaos_perturb()
        return self.gbest_pos, self.gbest_val, convergence

# ================== 遗传算法 ==================
class GA(Optimizer):
    """遗传算法"""
    def __init__(self, obj_func, dim, bounds, pop_size=50, max_iter=500,
                 crossover_rate=0.8, mutation_rate=0.1):
        super().__init__(obj_func, dim, bounds, pop_size, max_iter)
        self.cr = crossover_rate
        self.mr = mutation_rate
        
    def _select(self, fitness):
        """轮盘赌选择"""
        prob = 1 / (fitness + 1e-8)
        prob /= prob.sum()
        return self.pop_pos[np.random.choice(
            self.pop_size, self.pop_size, p=prob)]
    
    def _crossover(self, parents):
        """算术交叉"""
        children = parents.copy()
        for i in range(0, self.pop_size, 2):
            if np.random.rand() < self.cr:
                alpha = np.random.rand()
                children[i] = alpha*parents[i] + (1-alpha)*parents[i+1]
                children[i+1] = alpha*parents[i+1] + (1-alpha)*parents[i]
        return children
    
    def _mutate(self, children):
        """高斯变异"""
        mask = np.random.rand(*children.shape) < self.mr
        noise = np.random.normal(0, 0.1*self.range, children.shape)
        return np.clip(children + mask*noise, self.bounds[:,0], self.bounds[:,1])
    
    def optimize(self):
        self._init_population()
        convergence = []
        
        for _ in range(self.max_iter):
            # 评估
            fitness = np.array([self.obj_func(x) for x in self.pop_pos])
            
            # 选择
            parents = self._select(fitness)
            
            # 交叉
            children = self._crossover(parents)
            
            # 变异
            new_pop = self._mutate(children)
            
            # 更新
            self.pop_pos = new_pop
            current_best = np.min(fitness)
            if current_best < self.gbest_val:
                self.gbest_val = current_best
                self.gbest_pos = self.pop_pos[np.argmin(fitness)]
            
            convergence.append(self.gbest_val)
        
        return self.gbest_pos, self.gbest_val, convergence

# ================== 差分进化算法 ==================
class DE(Optimizer):
    """差分进化算法"""
    def __init__(self, obj_func, dim, bounds, pop_size=50, max_iter=500,
                 crossover_prob=0.8, scaling_factor=0.5):
        super().__init__(obj_func, dim, bounds, pop_size, max_iter)
        self.cr = crossover_prob
        self.sf = scaling_factor
        
    def _mutate(self):
        """差分变异"""
        mutants = np.zeros_like(self.pop_pos)
        for i in range(self.pop_size):
            idxs = [idx for idx in range(self.pop_size) if idx != i]
            a, b, c = self.pop_pos[np.random.choice(idxs, 3, replace=False)]
            mutants[i] = a + self.sf * (b - c)
        return mutants
    
    def _crossover(self, mutants):
        """二项式交叉"""
        cross_points = np.random.rand(*self.pop_pos.shape) < self.cr
        trial = np.where(cross_points, mutants, self.pop_pos)
        return trial
    
    def optimize(self):
        self._init_population()
        convergence = []
        
        for _ in range(self.max_iter):
            # 变异
            mutants = self._mutate()
            
            # 交叉
            trial = self._crossover(mutants)
            
            # 评估
            trial_val = np.array([self.obj_func(x) for x in trial])
            pop_val = np.array([self.obj_func(x) for x in self.pop_pos])
            
            # 更新
            self.pop_pos[trial_val < pop_val] = trial[trial_val < pop_val]
            current_best = np.min(pop_val)
            if current_best < self.gbest_val:
                self.gbest_val = current_best
                self.gbest_pos = self.pop_pos[np.argmin(pop_val)]
            
            convergence.append(self.gbest_val)
        
        return self.gbest_pos, self.gbest_val, convergence

# ================== 更复杂的多峰测试函数 ==================
def complex_rastrigin(x):
    """更复杂的Rastrigin函数"""
    return 10*len(x) + np.sum(x**2 - 10*np.cos(2*np.pi*x)) + np.sum(0.5*x**4)

def complex_griewank(x):
    """更复杂的Griewank函数"""
    term1 = np.sum(x**2 / 4000)
    term2 = np.prod(np.cos(x / np.sqrt(np.arange(1, len(x)+1))))
    return term1 - term2 + 1 + np.sum(0.1*x**4)

def complex_schwefel(x):
    """更复杂的Schwefel函数"""
    return np.sum(-x * np.sin(np.sqrt(np.abs(x)))) + np.sum(0.001*x**4)

def shifted_sphere(x):
    """高维Shifted Sphere函数"""
    shift = np.array([1.0, -2.3, 3.4, -4.5, 5.6, -6.7, 7.8, -8.9, 9.0, -10.1])
    return np.sum((x - shift[:len(x)])**2)

def noisy_ackley(x):
    """带噪声的Ackley函数"""
    noise = np.random.normal(0, 0.1, x.shape)
    return -20 * np.exp(-0.2 * np.sqrt(np.mean(x**2))) - np.exp(np.mean(np.cos(2*np.pi*(x + noise)))) + 20 + np.exp(1)

# ================== 实验执行 ==================
def main():
    # 测试函数
    test_functions = {
        'Complex Rastrigin': complex_rastrigin,
        'Complex Griewank': complex_griewank,
        'Complex Schwefel': complex_schwefel,
        'Shifted Sphere': shifted_sphere,
        'Noisy Ackley': noisy_ackley
    }
    
    # 算法配置
    algorithms = {
        'PSO': PSO,
        'CPSO': CPSO,
        'GA': GA,
        'ECPSO': ECPSO,
        'DE': DE
    }
    
    print("=== 算法对比测试 ===")
    results = {}
    
    # 基准函数测试
    for algo_name in algorithms:
        print(f"正在运行 {algo_name}...")
        optimizer = algorithms[algo_name](complex_rastrigin, dim=10, 
                                        bounds=[(-5.12,5.12)]*10)
        *_, curve = optimizer.optimize()
        results[algo_name] = curve
    
    # 收敛曲线可视化
    plt.figure(figsize=(12,6))
    for name, curve in results.items():
        plt.semilogy(curve, label=name, linewidth=2)
    
    plt.title("算法收敛曲线对比 (复杂Rastrigin函数)", fontproperties=zh_font, fontsize=14)
    plt.xlabel("迭代次数", fontproperties=zh_font)
    plt.ylabel("适应度值（对数尺度）", fontproperties=zh_font)
    plt.legend(prop=zh_font)
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.tight_layout()
    
    # 动态环境测试
    env = DynamicPeaks(dim=10, speed=0.1)
    dynamic_results = {}
    
    print("\n=== 动态环境测试 ===")
    for algo_name in algorithms:
        print(f"正在测试 {algo_name}...")
        tracker = algorithms[algo_name](env, dim=10, bounds=[(-5,5)]*10)
        errors = []
        for epoch in range(20):
            tracker.optimize()
            errors.append(np.linalg.norm(tracker.gbest_pos - env.peaks[env.current]))
            env.shift()
        dynamic_results[algo_name] = np.mean(errors)
    
    # 柱状图可视化
    plt.figure(figsize=(10,6))
    names = list(dynamic_results.keys())
    values = list(dynamic_results.values())
    plt.bar(names, values, color=['#4C72B0','#DD8452','#55A868','#C44E52','#8172B2'])
    
    plt.title("动态环境平均追踪误差", fontproperties=zh_font, fontsize=14)
    plt.ylabel("平均位置误差", fontproperties=zh_font)
    plt.xticks(fontproperties=zh_font)
    plt.grid(True, axis='y', linestyle='--', alpha=0.7)
    
    # 多峰函数测试成功率
    success_rate = {}
    for algo_name in algorithms:
        print(f"\n=== 测试 {algo_name} 在复杂多峰函数上的成功率 ===")
        successes = 0
        for _ in range(10):
            optimizer = algorithms[algo_name](complex_rastrigin, dim=10, 
                                           bounds=[(-5.12,5.12)]*10)
            _, best_val, _ = optimizer.optimize()
            if best_val < 1e-3:
                successes += 1
        success_rate[algo_name] = successes / 10
    
    # 成功率柱状图
    plt.figure(figsize=(10,6))
    names = list(success_rate.keys())
    values = list(success_rate.values())
    plt.bar(names, values, color=['#4C72B0','#DD8452','#55A868','#C44E52','#8172B2'])
    
    plt.title("复杂多峰函数全局寻优成功率", fontproperties=zh_font, fontsize=14)
    plt.ylabel("成功率", fontproperties=zh_font)
    plt.xticks(fontproperties=zh_font)
    plt.ylim(0, 1)
    plt.grid(True, axis='y', linestyle='--', alpha=0.7)
    
    # 参数敏感性分析
    print("\n=== 参数敏感性分析 ===")
    param_results = {}
    param_values = [0.5, 0.729, 0.9]
    for omega in param_values:
        optimizer = ECPSO(complex_rastrigin, dim=10, bounds=[(-5.12,5.12)]*10)
        optimizer.omega = omega
        *_, curve = optimizer.optimize()
        param_results[f'omega={omega}'] = curve
    
    # 参数敏感性曲线
    plt.figure(figsize=(12,6))
    for name, curve in param_results.items():
        plt.semilogy(curve, label=name, linewidth=2)
    
    plt.title("参数敏感性分析 (omega)", fontproperties=zh_font, fontsize=14)
    plt.xlabel("迭代次数", fontproperties=zh_font)
    plt.ylabel("适应度值（对数尺度）", fontproperties=zh_font)
    plt.legend(prop=zh_font)
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.tight_layout()
    
    # 高维测试
    print("\n=== 高维测试 ===")
    high_dim_results = {}
    for algo_name in algorithms:
        print(f"正在测试 {algo_name} 在20维复杂Rastrigin函数上...")
        optimizer = algorithms[algo_name](complex_rastrigin, dim=20, 
                                        bounds=[(-5.12,5.12)]*20)
        *_, curve = optimizer.optimize()
        high_dim_results[algo_name] = curve
    
    # 高维测试收敛曲线
    plt.figure(figsize=(12,6))
    for name, curve in high_dim_results.items():
        plt.semilogy(curve, label=name, linewidth=2)
    
    plt.title("20维复杂Rastrigin函数收敛曲线", fontproperties=zh_font, fontsize=14)
    plt.xlabel("迭代次数", fontproperties=zh_font)
    plt.ylabel("适应度值（对数尺度）", fontproperties=zh_font)
    plt.legend(prop=zh_font)
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.tight_layout()
    
    # 噪声环境测试
    print("\n=== 噪声环境测试 ===")
    noisy_results = {}
    for algo_name in algorithms:
        print(f"正在测试 {algo_name} 在带噪声的Ackley函数上...")
        optimizer = algorithms[algo_name](noisy_ackley, dim=10, 
                                        bounds=[(-32.768,32.768)]*10)
        *_, curve = optimizer.optimize()
        noisy_results[algo_name] = curve
    
    # 噪声环境收敛曲线
    plt.figure(figsize=(12,6))
    for name, curve in noisy_results.items():
        plt.semilogy(curve, label=name, linewidth=2)
    
    plt.title("带噪声的Ackley函数收敛曲线", fontproperties=zh_font, fontsize=14)
    plt.xlabel("迭代次数", fontproperties=zh_font)
    plt.ylabel("适应度值（对数尺度）", fontproperties=zh_font)
    plt.legend(prop=zh_font)
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.tight_layout()
    
    plt.show()

if __name__ == "__main__":
    main()